function [outFEM] = wrapper_fem (  x_nodes , x_samples , conectivity, gPts) %, order

sPts = size(x_samples,1);
nElements = size (conectivity,1);

p_a  = {zeros(sPts,1)};
dp_a = {zeros(sPts,1)};

ig=0;

for elem= 1:nElements
    Xa    = x_nodes(conectivity(elem,:),:); %vertices
    for gp =1 : gPts
        ig=ig+1;
        [p_a{ig} dp_a{ig}] = TriangleBasisFunctions (Xa , x_samples (ig,:) );
    end
end

outFEM.p_samp = p_a; 
outFEM.dp_samp = dp_a;


end
